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ABSTRACT 

We study the stability of standing shock waves in advection-dominated accre- 
tion flows into a Schwarzschild black hole by 2D general relativistic hydrodynamic 
simulations as well as linear analysis in the equatorial plane. We demonstrate 
that the accretion shock is stable against axisymmetric perturbations but be- 
comes unstable to non-axisymmetric perturbations. The results of dynamical 
simulations show good agreement with linear analysis on the stability, oscillation 
and growing time scales. The comparison of different wave-travel times with the 
growth time scales of the instability suggests that the instability is likely to be 
of the Papaloizou-Pringle type, induced by the repeated propagations of acous- 
tic waves. However, the wavelengths of perturbations are too long to clearly 
define the reflection point. By analyzing the non-linear phase in the dynamical 
simulations, it is shown that quadratic mode couplings precede the non-linear 
saturation. It is also found that not only short-term random fluctuations by 
turbulent motions but also quasi periodic oscillations take place on longer time 
scales in the non-linear phase. We give some possible implications of the insta- 
bility for quasi periodic oscillations (QPOs) and the central engine for gamma 
ray bursts (GRBs). 

Subject headings: GRB , central engine , SASI , QPO , accretion disk , GRHD 

1. Introduction 

Accretion flows imbedding a shock wave have attracted much attention of researchers. 
Hydrodynamic instabilities of shocked accretion flows may explain the time variability of 
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emissions from many black hole candidates, since a shock wave is a promising mechanism of 
transforming the gravitational energy into radiat ion. The poss i bility t hat a shock exists in 
black hole accretion disks was first sugg ested by lliawlev et all Jl984al lbl) even though they 
did not make clear the essential condition for its existence. The possible structure of the 
shocked accretion flows in equatorial plane was described by Fukuel (jl987l ). who suggested 
the importance of the multiplicity of critical points for the existence of the standing shock. 
Multiple critical points exist only for appropriate values of the injection parameters such as 
the specific angular momentum and Bernoulli constant. Note also that there are generally 
two possible shock locations, which are referred to as the inner and outer shocks. 

The stability of the standing shock wave in the accr etion flows has been also investigated 
by many authors both analytically and numerically. iNakayamal (119941 119951 ) showed by 
linear analysis in the equatorial plane that if the post-shock matter are accelerated, the 
flow is unstable against radial perturbations, which is true for both Newtonian or general 
relativistic dynamics. As a result of this theorem, we know for the black hole accretion 
that the inner shock becomes generally unstable against radial perturbations and that non- 
rotating steady accretion flows to a black hole cannot have a stable standing shock wave in it. 
These feat ures were also observed in ID axisymmetric simulatio ns for a pseudo-Newtonian 
potential Jchakrabarti fc Molten] Il993l : iNobuta fc Hanawalll994l ). 



Recently, iFoglizzol (120011 . |2002| ) pointed out by linear analysis that the outer shock wave, 
which is stable against radial perturbations, is in fact unstable to non-radial axisymmetric 
perturbations. He argued that the advective-acoustic cycle could be responsible for the in- 
stability. In this mechanism the velocity and entropy fluctuations initially generated at the 
shock are advected inwards, producing pressure perturbations, which then propagate out- 
wards, reach the shock and generate entropy and velocity fluctuations there, thus repeating 
the cycle with an increased amplitude. The same instability appears to be working in the 
accretion flows onto a nascent proto neutron star in the super nova core, in which l arge scale 
oscillations of the shock wave of £ = 1 nature are observed (IBlondin et al.ll2003l ). where I 
stands for the index in the spherical harmonic functions Y^. Although the mechanism is still 
contr oversial (jOhnishi et al.ll2005l : iFoglizzo et al.ll2006l ; IBlondin fc Mezzacappall2006l ; lLaming 
20071 ). the so-called standing accretion shock instability or SASI is currently attracting much 
attention as a promising explanation of asymmetr ic explosion of supernova as well as young 
pulsars' proper mot ions (ISchek et al.l 120041 . 120061 ) and spins (IBlondin fc Mezzacappal 120071 ; 
Iwakami et al.l 120081 ) . 



As for the stability of the shock wave against non-axisymmetric perturbations, iMolteni et al 



(119991 ) did 2D simulations of an adiabatic shocked accretion flow by using the pseudo- 
Newtonian potential and found a non-axisymmetric instability. They showed that the shock 
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instability is saturated at a low level and a new quasi-ste ady asymmetri c confi g uration is 



realize d. To investigate the mechanism of this instability, iGu &: Foglizzol (120031 ); iGu &: Lu 



( 120051 ) performed linear analysis for both isothermal and adiabatic flows. They concluded 
that the instability seems to be of the Papaloizou-Pringle type and the repeated propaga- 
tions of acoustic waves between the corotation radius and the shock surface are a driving 
mechanism. These conclusions are based on the WKB approximation and the comparison 
between the growth rate and the acoustic cycle period. 



It is also noted that lYamasaki fc Foglizzol (120071 ) investigated the linear stability of the 



shocked accretion flows onto the proto neutron star against non-axisymmetric perturbations 
in the equatorial plane. They demonstrated that the counter-rotating spiral modes are sig- 
nificantly damped, whereas the growth rate of the corotating modes is increased by rotation. 
They also claimed that the instability is not of the Papaloizou-Pringle type, since the stability 
is not affected by the presence or absence of the corotation point. Instead, they suggested 
the advective-acoustic cycle based on the WKB analysis. The purely acoustic cycles was 
found to be stable. 

As mentioned above, although many efforts have been made for clarifying the non-radial 
instability, the complete understanding of the mechanism is still elusive. It is noted that 
the unperturbed accretion flows should be treated properly, since they strongly affect the 
instability. In black hole accretions, the gravitation is one of the main factors to determine 
the flow features such as sonic points and shock locations. So far, however, the shock stability 
in the accretion flows to black holes has been investigated under the Newtonian or pseudo- 
Newtonian approximation and there has been no fully GR treatment. As shown later, the 
range of injection parameters that allows the existence of a standing shock wave is changed 
substantially when GR is fully taken into account. 

In this paper, we investigate fully general relativistically the stability of the shock in 
the advection-dominated accretion flows into a Schwarzschild black hole by using both linear 
analysis and non-linear dynamical simulations. In so doing, we consider only the equatorial 
plane, assuming that the ^-component of four-velocity (u e ) and all the ^-derivatives are 
vanishing. We use the spherical coordinates (r, 6, <p) in the following. The evolution of the 
metric is not taken into account, which is justified if the mass of accretion flow is much 
smaller than the black hole mass. We show that the shock is indeed unstable against non- 
axisymmetric perturbations and a spiral arm structure is formed as the instability grows. 
We discuss the instability mechanisms comparing various time scales. Finally, we mention 
possible implications of our findings for gamma-ray bursts (GRBs) and black hole quasi- 
periodic oscillations (QPOs). 

This paper is organized as follows. In section 2, we describe the steady axisymmetric 
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accretion flows with a shock in them. In section 3, we present the formulation of linear 
analysis. The numerical method for dynamical simulations is explained in section 4. The 
main numerical results and their analyses are given in section 5. The implication for GRBs 
and Black Hole QPOs are mentioned in section 6. Finally we summarize and conclude the 
paper in section 7. 



One of the key features of the accretion flows to black holes is that the inflow velocity 
is supersonic at the event horizon. This immediately means that there should be at least 
two sonic points if a steady shock wave exists in the accretion flows, since both the pre- and 
post-shock flows are transonic. This is in sharp contrast to the accretions onto a neutron 
star, in which the post shock flow is subsonic. One of the consequences of this theorem 
is the fact that spherical adiabatic accretions into Schwarzschild black holes are unable to 
have a steady shock wave in them, since they have a single sonic point. As long as rotating 
accretion flows in the equatorial plane are concerned, the locations of these sonic points are 
determined by the adiabatic index and injection parameters, such as the Bernoulli constant 
and specific angular momentum. Note that the accretion rate is irrelevant for the locations 
of sonic points and standing shocks. 

The basic equations are the relativistic continuity equation and equation of energy- 
momentum conservations: 



where the Greek indices represent the spacetime components. As already mentioned, we con- 
sider the accretions only in the equatorial plane, assuming ^-component of velocity (u e ) and 
all 9 derivatives are vanishing. Then the basic equations are reduced to ordinary differential 
equations with respect to the radial coordinate: 



2. Axisymmetric Steady Accretion Flows with a Shock 



2.1. The multiplicity of sonic points 



(Pou% = 0, 
(T»% = 0, 



(1) 
(2) 
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Note t h at these treatme nts are slightly different from those by lMolteni et al.l (119991 ) ; iGu fc Foglizzo 
( 120031 ) ; IGu fc Lul (120051 ) . They employed the cylindrical coordinates and integrated out the 
vertical structures, thus considering the accretion flows only in the equatorial plane. On the 
other hand, we use the spherical coordinates in this paper, since it is mathematically more 
convenient for the fully general relativistic approach. Because of this difference, however, 
the obtained accretion flows would be different from the ones considered in this paper if the 
general relativity were taken into account in their formulations. 

The left panel of Figure Q] shows the locations of the sonic points in Schwarzschild 
blac k hole as a function of the Bernoulli constant and specific angular momentum (see 
also [lu| ( 1986 )). As shown in the figure, there are indeed two or three sonic points for 
some combinations of the Bernoulli constant and angular momentum. The maximum and 
minimum specific angular momenta are A 

max ~ 4.0iW* and A m j n ~ 3.2ilf*, respectively, for 
the Bernoulli constant E = 1.004, for example. It should be noted that one of the most 
important differences between the full GR and the pseudo-Newtonian treatments is the 
range of the injection parameters that allows the existence of multiple sonic points. As the 
Bernoulli constant becomes smaller, the maximum specific angular momentum gets larger 
without limit for the pseudo Newtonian case whereas it is bounded for the full GR case. 
For the Bernoulli constant that may be typical for massive stellar collapse, e.g. E = 1.003, 
the maximum specific angular momentum is larger by about 60% for the pseudo Newtonian 
approximation than for the full GR treatment. The reason for this difference is that the 
gravity nearby the black hole is too strong in the pseudo-Newtonian approximation. It 
should be pointed out that the region of the injection parameters that allows multiple sonic 
points is not very wide and the above-mentioned difference may be important in considering 
the implications for astrophysical phenomena. 



2.2. The locations of standing shock waves 

The sonic points discussed in the previous subsection correspond to the so-called critical 
points in dynamical system. It is known that the innermost and outermost critical points 
are of the saddle-type, while the middle critical point is of the center-type. Hence, the 
transonic accretion flows can be constructed only for the former two critical points. These 
two transonic flows have the same Bernoulli constant and angular momentum but different 
entropies, and they can be connected by a standing shock wave, where the Rankine-Hugoniot 
relations hold: 



[p u^]^ = 0, 



(7) 
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[T^\l v = 0, (8) 

where is a 4- dimensional vector normal to the shocked surface, and is set to be = 
(0, 1, 0, 0) in axisymmetric steady flows (see also Eq. (123]) for the non-axisymmetric per- 
turbations). We use the notation of [Q] = Q + — Q_, where the subscript (+) represents 
a post-shock quantity and (— ) a pre-shock quantity. The Bernoulli constant and specific 
angular momentum are defined respectively as 

E = -hu t , (9) 

A EE (10) 

Ut 

Since the Bernoulli constant, specific angular momentum and mass flux are unchanged across 
the shock, we only need to consider the radial component of energy momentum conservation 
across the shock. Then the Rankine-Hugoniot relation in axisymmetric steady flows can be 
written as: 

[p hu r u r +p} = 0. (11) 

In general, there are two possible shock locations, to which we refer as the inner and 
outer shocks. This is apparent in the right panel of the Figures [TJ in which we show the 
parameter region that allows multiple shock locations. It is w ell known, howev er, that the 



inner shock is unstable against axisymmetric perturbations ( iNakayamal Il995l ). which we 
have confirmed by our own linear analysis. We have also done numerical simulations with a 
single grid point in the azimuthal direction, thus suppressing non-axisymmetric modes, and 
observed that the inner shock is either swallowed by a black hole or moved outwards to be 
the outer shock after radial perturbations are imposed. On the other hand, we have seen 
that the outer shock is stable against radial perturbations even if the perturbation amplitude 
is not necessarily small. In the following, we will consider only the outer shock. 

We have constructed several axisymmetric steady accretion flows with an outer shock for 
different combinations of the adiabatic index and injection parameters, which are summarized 
in Table [TJ 



3. Linear Analysis of non-axisymmetric shock instability 

Here we give the basic equations and boundary conditions for the linear analysis. The 
obtained eigen values are later compared with the numerical simulations, and the eigenstates 
are employed to impose the initial perturbations. 
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The basic equations are the linearized relativistic continuity and energy-momentum 
tensor conservation equations (see Eqs. ([I]) and (j2J)). We again neglect the ^-component 
of velocity and all the ^-derivatives and consider the equatorial plane only Under this 
assumption, the system equations are written as follows: 

drf = po(0) V(o) { Po(i)u m ° + po(o) {^ t(1) ~ mu^) }, (12) 

9 r q = zJrr m ( u Pi + Po(o)U m h {0 )u m (rq) , (13) 

u t(0) 

9rV {1) = (14) 

u t{0) 

where S is entropy and the following notations are used: 

_ po(i) « r(1) , s 

; = P^ + ^)' (16) 

- ^ + ^, (17) 

tyo) Wi(0) 

V(i) = a; ( h^u^o) + fyo)iV(i) ) + m ( fyi)Ut(o) + fyo)«t(i) ) , (18) 
u 4>(o) 

Following the standard procedure of linear stability analysis, the perturbed quantities are as- 
sumed to be proportional to e - lw *+ im< A. All perturbed quantities are calculated from /, q, V(i) 
and S(i). V(i) and can be integrated analytically as 

/ u *(o) \ 

% = Vfo exp z— T o- , (20) 



5(i) = 



R - V M r(0) 
„t(0) 



R exp(z-^a), (21) 



where the subscript 'R' denotes the value evaluated at r = R. Thus, we need to integrate 
numerically only two Eqs. ffl2l) and (fT3l) . 

These linearized equations can be solved with appropriately setting boundary conditions, 
which are imposed at the shock surface and inner sonic point. In this study, we assume that 
the perturbations are confined in the post-shock region and the pre-shock region remains 
unperturbed. As a result, the outer boundary condition is set at the shock surface. We 
express the shock radius as follows: 

R s h = Rsh(o) + V exp ( -iut + inup ) , (22) 
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where rj denotes the initial amplitude of shock displacement. Defined as the 4-dimensional 
vector normal to the shock surface, l v can be given as 



Z„ = (iujrie- iult+im ' t ' J 1, 0, -imr]e- iujt+im ^) 



(23) 



Using these relations in the Rankine-Hugoniot relations expressed in general as [Q^ \ = 0, 
we can write Eqs. © and (|HJ) as 



[Q m }turi- \Q m ]imr l + [Q r(1) 
where the following notations are employed: 

Qrto) = Q m 

q M (i) = qh(i) +r) I "qm(o) 



0. 



HO) 



Rsh(O) 



d 



dr 



As mentioned above, we assume that the pre-shock quantities are unperturbed, i.e 
0. The explicit forms of these equations are given by 



(po(o)« r(1) ) + +(po(i)« r(0) ) + + A = 0, 



[p hu t u r ) { 



(o) / Po(i) . h(i) u 



t(i) 



Po(o) 



+ 



h 



+ 



(0) 



( Po h(u r ) 2 y 



2\(0) / Po(l) . h(i) 



u 

,r(l) 



t(0) 



+ 



r(l) 



r(0) 



B = 0, 



+ 



(pohu+vTy 
with the following definitions: 



Po(0) fyo) 

(o) / Po(i) , h(i) u 



u r(l) \ i 
+ 2^J+ P(1) ^} + + C = I 



0(i) 



P0(0) fyo) w 



0(0) 



+ 



u 



r(l) 



r(0) 



D = 0, 



.4 



d 
dr 



Po(o) 



r(0) 



d 
dr 



P0(0)M 



r(0) 



-im V \(po { o ) u^) + -(p oi o ) u^) 



Ho)" 



B = 



d 
dr 



p (o)W (0) « r(0) 



dr 



p m h u^u^) U 



Po(o)^(o) K (0) ) 2 +P (0) ^ 



C = 



(o) 



d_ 

dr 



d 
dr 



(24) 

(25) 
(26) 



(27) 
(28) 

(29) 

(30) 



(31) 



(32) 



P0(0)^(0)K (0) ) 2 +pW^)) }r? 
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+ioor] | (po^V)| 0) - (p /^V) W } 

-imr] | {poh^u^f - { Po hu r u*) m } , (33) 
D - {(^Po(o)^ ( %^) + -(^o ( o)^%^) 
+z^{(p /mV)j 0) - (p /mV) W } 

-z'mtj { (p /i K) 2 +p^) (0) - (pah (u*) 2 +pg**)®}. (34) 

As for the inner boundary condition, on the other hand, a regularity condition is imposed 
at the sonic point. By combining Eqs. ffl~2l) to (fT9l) . we obtain the differential equation for 
u r , which is written generally as 

F(u r ) r = G. (35) 
The linearized form of this equation becomes 

Km) r = G( " \ m "™' , (36) 

^(0) 

where the explicit forms of F(o), Fm and are 



F 



(o) = Po(o)h( )g rr {u r( -°h r{0 ) - (& s (o)) 2 (l + M r(0) M r(0 )) | (37) 
F {1) = p (i)/i(o) K (0) ) 2 + p (o)/i(i) K {0) ) 2 + 2p 0(0) /i (0) ?/ (( V (1) 

-r/ r {p (1) (1 + u r «V( )) + 2p (0)? /(V (0) }, (38) 

G(l) = r|V r ir + ^ {(P(l)M r ( ) +P(0)M r( l)) (1 +M r ( )M r(0) ) +2p (0 ) (M r(0 )) 2 « r(1) } 

+ ^{po(i)^ ( o)^ (0) + Po(o)/i(i)W r(0) + Po(o)^(o)^ (1) }{^ (n r(0) ) 2 + ^ K (0) ) 2 + <? tt , r (m' (0) ) 2 } 

+ PO(0)^(0)« r(0) {w« r(0) « r(1) + Q^r^U*^ + ^rtl^TI^} 
+ ^P0(0)^(0)W r(0) M* (0) Mr.(l)O- 

+ iTp (0) (1 + u r[0) u r( - 0) ) (-um*« + m^W) 

-iu m p (1) a. (39) 

6 s (o) and T are the unperturbed sound velocity in the comoving frame and the adiabatic index, 
respectively. Since F( ) vanishes at the sonic point, we impose the regularity condition there: 

G ( i) - F (1) u r(0) , r = 0. (40) 
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4. Numerical Method and Models 

The growth of the initial perturbations is computed with a multi-dimensional general 
relativistic hydrodynamics code, which is based on a modern technique, the so-called high- 



resolution central scheme (IKurganov fc Tadmorll2000l ). The details of the numerical scheme 
are given in Appendix. 

We use the Kerr-Schild coordinates with the Kerr parameter being set to be zero, since 
they have no coordinate singularity at the event horizon and we can put the inner boundary 
inside the event horizon. This is very advantageous for numerical simulations. We employ 
a T-law EOS, p = (T — 1) p e, where p and p and e are the pressure, rest-mass density, and 
specific internal energy, respectively. 

The computational domain is a part of the equatorial plane with 1.5M* < r < 200M* 
and the computation times are ~ 6 x 10 4 M^, where M* denotes the black hole mass and 
the unit with c = G = 1 is used. We employ 600(r) x 60(0) grid points. The radial grid 
width is non-uniform with the grid being smallest, (Ar = 0.1M*), at the inner boundary 
and increasing geometrically by 0.34% per zone toward the outer boundary. 

The initial perturbation modes are summarized in Table [TJ For Models Ml to M12 
we add the m = 1 mode, where m stands for the azimuthal mode number in e m< ^. For 
Models Mlm2 and Mlm3, on the other hand, the m = 2, 3 modes are initially imposed, 
respectively. These models are meant for the study of the initial-mode dependence. In order 
to investigate the initial-amplitude dependence, we also run Models MlalO and MlalOO, 
whose initial amplitudes are 10% and 100%, respectively. The unperturbed flows for Models 
Mlm2, Mlm3, MlalO and MlalOO are common to those of the other models. 

The radial distributions of these modes are obtained by the linear analysis. This is 
important to compare the linear growth of purely single modes between t he lin ear analysis 



and numerical simulations, and is one of the differences from lMolteni et al.l (119991 ) . We choose 
the most unstable mode for each azimuthal wave number except for MlalOO, in which the 
following perturbation is employed to avoid a velocity larger than the light velocity: 



^ ta {l + sin(0)}, (41) 



where v r sta is the unperturbed radial velocity. The initial perturbations are added to the 
whole post-shock region, except for Model MlalOO, in which the perturbation is imposed to 
the region between the inner sonic point and the shock again to prevent the flow velocity 
from being larger than the light velocity. 

In the following, we set the black hole mass to be M* = 3M , where M & is the solar 
mass. We have in mind here some applications to astrophysical phenomena such as GRBs 
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and QPOs. Since the formulation is dimensionless, the scaling is quite obvious. 



5. Numerical Results by GRHD 



In this section, we describe the time evolutions of non-axisymmetric instability obtained 
by fully dynamical GRHD simulations. In the following analysis, we frequently employ the 
mode decomposition of the shock surface by the following Fourier transform: 



Jo 

where R s h (0, t) and a m (t) are the radius of the shock wave as a function of (ft and t and the 
amplitude of mode m as a function of t, respectively. 



Figures [2] and [3] show the temporal evolutions of velocity and entropy for the baseline 
model Ml, respectively. The perturbation grows exponentially at the beginning and the 
shock wave is deformed according to the imposed m — 1 mode, for which the deformed 
shock surface rotates progressively, that is, in the direction of the unperturbed flows. Then 
the shock radius, or the m = mode, starts to grow. After several revolutions, a spiral 
arm develops and the instability is saturated with more complex structures. In this non- 
linear regime, several shocks are generated and collide with each other. As a result of these 
interactions, the original shock oscillates radially. 

As given in Table [2} the saturation levels of shock radius, or the m = mode, differ 
widely among models. For example, the left panel of Figure H] shows the time evolution of 
m = mode for Model Ml. Both large- and small-amplitude oscillations can be seen, which 
correspond to the periods of ~ 100 ms and ~ 20 ms, respectively. These two kinds of ax- 
isymmetric oscillations are also found in other models. Their characteristics are summarized 
in Table [2J The growth of m = mode in the non-linear phase is strongly dependent on the 
mach number (see Tables [1] and [2] ) . i.e., a stronger shock tends to become more unstable 
against non-axisymmetric perturbations. Although there is no entry to the maximum am- 
plitude of m=0 mode in Table 2 for Model M5, this is not an exception. In fact, the shock 
is so unstable in this model that it leaves the computational domain soon after the addition 
of non-axisymmetric perturbations. 

On the other hand, the oscillation periods of the non-axisymmetric modes are much 
shorter in general than those of the axisymmetric ones. Indeed, the right panel of Figure H] 




(42) 



5.1. 



Basic Features 
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shows the time evolutions of m = 1,2,3 modes for Model Ml. The typical periods range 
from several to a few dozen milli-seconds. 

In the non-linear phase, the dominant mode is m — 1 for almost all models. Remarkably, 
although the initially added perturbations are not the m = 1 mode in Models Mlm2 and 
Mlm3, the non-linear mode couplings lead eventually to the dominance of the m = 1 mode 
(see the upper panels of Figure [9]). 



5.2. Comparison with linear analysis 

Equations (TP21) to (fl5|) with the boundary conditions at the shock surface and sonic point 
are solved numerically to find eigen modes. Figure [5] shows the real and imaginary parts of 
eigen frequencies for some of the m = 1,2,3 modes for Model Ml. They are all unstable 
non-axisymmetric modes. We find, on the other hand, that the axisymmetric perturbations 
are st able, which is con sistent with both the present dynamical simulations and the previous 



work (jNakayamal Il995l ) . The oscillation periods, which correspond to u r , are 1.5 ~ 3.7 ms 
for the most unstable mode in each m-sequence for Model Ml, whereas the growth times, 
which are obtained from uii, are 2.6 ~ 3.2 ms. 

Figure [6] shows the comparison of the amplitudes obtained by the linear analysis and 
dynamical simulations. As can be clearly seen, the oscillation period and growth time are in 
good agreement between them for the first 10 ms, which is the linear growth phase. It is also 
found that after 10 ms, the result of the dynamical simulation starts to deviate from that of 
the linear analysis, which indicates the beginning of the non-linear phase. The amplitude is 
saturated and the oscillation period gets slightly longer. According to Figure El which shows 
the evolution of axisymmetric m = mode, we find that it starts to grow at t ~ 10 ms, 
leading to the increase of the average shock radius. This is the reason why the oscillation 
period becomes longer in the non-linear phase. 

As mentioned in the previous section and summarized in Table |2J the dominant modes 
in the non-linear phase are progressive, that is, the deformation pattern rotates in the same 
direction as the unperturbed flow. This is also true for all the linearly unstable modes. 
In fact, the linear analysis shows that there is no unstable mode for m < 0. Note that 



Yamasaki fc Foglizzol (120071 ) demonstrated by linear analysis for the accretion onto a neu- 
tron star that the progressive modes are enhanced and retrogressive ones are suppressed by 
rotation. 



As the specific angular momentum of the unperturbed flow becomes larger, the distance 
between the shock wave and the inner sonic point gets greater (See Models Ml and M4 in 
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Tabled]). According to the linear analysis, the number of unstable modes also increases, 
whereas the growth time of the most unstable modes becomes longer. This suggests that 
the larger distance tends to stabilize the non-axisymmetric instability though it is not the 
only factor for the shock instability. In fact, it is also found that the growth rate of unstable 
modes is affected by the shock strength, that is, stronger shocks tend to be more unstable. 

Finally we show in Figure [8] the most unstable modes for different m's in Model Ml. As 
is clear, the most unstable of all is the m = 4 mode and the modes with m > 12 are stable 
for this model. It is also interesting to note that the real part of eigen frequency for the 
most unstable modes becomes larger as the mode number m is higher (see the left panel of 
Figure [Sj), while the pattern frequency u r /m becomes smaller (the right panel of Figure [8]). 

5.3. Dependence on the initial perturbations 

By comparing Models Ml, MlalO and MlalOO, we find that the qualitative feature 
of dynamics are almost the same. In the linear phase, the growth of mode amplitudes is 
unaffected by the initial condition. Only the duration of the linear phase become shorter 
as the initial amplitude is larger as expected. The saturation levels do not differ very much 
among three models. In fact, the non-linear phase is rather chaotic and forgets the difference 
in the initial condition. It is also important to point out that in spite of the large initial 
amplitude for Model MlalOO, the shock continues to exist. This implies that if the injection 
parameters are appropriate, a standing shock will exist oscillating violently in the accretion 
flows into black holes. 

Next we show what happens if we impose initially the m = 2 or 3 mode instead of the 
m — 1 mode, comparing Models Ml, Mlm2 and Mlm3. In the upper panels of Figure [91 we 
show the time evolutions of the amplitudes for various modes in Models Mlm2 (left panel) 
and Mlm3 (right panel). 

We first pay attention to the evolution up to ~ 150 ms, where the difference is most 
evident. Note that the linear phase lasts only for ~ 10 ms and the non-linear phase thereafter 
is the focus here. In the left panel, we see the growth and saturation of m = 4 mode in 
addition to the original m = 2 mode. On the other hand, the m = 6 mode is formed to 
grow to the saturation in the right panel. Note that the m = mode is also generated 
in these models and will be discussed later. These models are produced by the non-linear 
mode-couplings, which are of quadratic nature, and the other modes with odd m for Model 
Mlm2, for example, are not generated. 

After ~ 150 ms, however, other modes also emerge and grow to be saturated. These 
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modes are probably generated by numerical noises, having much smaller amplitudes initially 
and spending longer time in the linear phase. After the saturation of the modes, the dynamics 
is almost identical for the three models. 

In the lower panel of Figure [9j we show the time evolution of the amplitudes of m = 
mode for the three models. As mentioned above, the modes are produced by the quadratic 
mode couplings of the initially imposed modes. Up to ~ 150 ms, the evolutions are quite 
different among them. For Model Ml, the amplitude grows to be 4 times larger than the 
initial value in ~ 50 ms and oscillates violently thereafter. On the other hand, the maximum 
amplitudes are much smaller for Models Mlm2 and Mlm3, and the following oscillations 
have much smaller amplitudes. In fact, for Model Mlm3 the amplitude becomes almost 
constant after ~ 50 ms. The saturation level is much lower than that of Model Mlm2. It is 
interesting to point out that even though the m = 3 mode is the unstable in the linear phase 
(see Figure EJ), the average shock radius, or the m = mode, is most strongly affected by 
the m = 1 mode. After ~ 150 ms, all modes are saturated and the behavior of the m = 
mode becomes almost identical among the different models. 



5.4. Dependence on the adiabatic index 

In this paper, we employ the simple T-law EOS and so far we have discussed only the 
case with T = 4/3. In reality the EOS will not be so simple and the adiabatic index may not 
be constant. In order to infer the differences that the EOS may make, we vary the adiabatic 
index in the the T-law EOS and see the changes in this subsection. 

It is the unperturbed accretion flows that are most affected by the change of adiabatic 
index. It is found that as the adiabatic index becomes larger, both the specific angular 
momentum and Bernoulli constant that allow the existence of a standing shock wave gets 
smaller. This is understood as follows. The structure of unperturbed accretion flows and 
hence the existence of shock are determined by the balance between the attractive gravity 
and the repulsive centrifugal force and pressure. As the EOS becomes harder, the pressure 
gets larger and, as a result, the centrifugal force can be reduced. For the same specific 
angular momentum, on the other hand, the Bernoulli constant can be smaller for harder 
EOS's. Note that the Bernoulli constant is a measure of the matter temperature at infinity. 

The instability itself is also affected by the change of adiabatic index, since it depends on 
the structure of unperturbed accretion flows. The results are summarized in Tables [2] and [3j 
Although it is difficult to find a systematic trend here, the saturation amplitude of the m = 
mode appears to be correlated with the Mach number: as the Mach number becomes larger 
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by the change of adiabatic index, the saturation level gets enhanced. It is more important to 
understand here, however, that the instability does not change qualitatively in spite of the 
relatively large variation in the injection parameters that allow the existence of the standing 
shock wave. More thorough investigations of the EOS dependence will be a future task. 



5.5. Instability Mechanisms 



To gain some insight into the instability mechanism, we calculate some time scales as 
follows: 



T d-c(c) ■ 
1~c—c(c) 
T~d-c{s) '- 
1~c—c(s) 



r s h 



cor 



r cor 



1 11 S O 



1 1 

C7 + Cj 



dr. 



(43) 
(44) 
(45) 
(46) 



where r inso is the radius of inner sonic point and Cf are the outgoing (+) and ingoing (-) 
sound velocities, respectively, and are given in the observer's frame for the Schwarzschild 
geometry as, 



C 



(l-a M V±{{(l-6 s 2 ) M V} 2 -[(l-6 s 2 ) (u*) 2 - 6 s V*][(l-fo s 2 ) K) 2 - bfgT] 



(1-6 2 ) («')' - b*g« 



.(47) 



Here b s denotes the sound velocity in the comoving frame. The corotation point of the 
perturbation is defined as 



LO r - m— = 

u f (r cor ) 



(48) 



and its radius is expressed as y. This investigation is inspired by the previous works 



flGu fc Foglizzcl 12003 



Gu fc Lul 120051 ) mentioned in the introduction. These time scales and 



the radius of corotation point for all the models are summarized in Table El For comparison, 
we list the oscillation and growth time scales, which are obtained by linear analysis for the 
most unstable mode. 



Figure [10] compares the growth times with the cycle periods given by Eqs. (T43l) to (T46l) 
for all the models. It is found that the periods of acoustic-acoustic cycle are closer to 
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the growth times than those of the advective-acoust i c cycle, which appears to support the 



claim by the preceding papers (IGu fc Foglizzol 120031 ; iGu Lul 120051 ) that the instability is 



of the Papaloizou-Pringle type. It is important, however, to point out that we can not 
identify the reflection point clearly (see Figure [TDJ, in which the left (right) panel adopts 
the corotation (inner sonic) point as the inner reflection point). This is mainly because the 
wavelengths of the perturbations are rather long. In fact, we estimate the wavelength of 
acoustic perturbation by 

\ w = min (^^j , (49) 

which should be much smaller than the scale height of the unperturbed flows for the justi- 
fication of WKB approximations. The wavelength of the dominant unstable mode for each 
model is given in Table [3] According to this estimation, it is comparable or longer than 
the scale height. Hence the WKB approximation is not justified at least for these models. 
Indeed, the reflection point of waves lose its meaning. Incidentally, the WKB approximation 
may be applicable to higher harmonics. In fact, there are sequences of unstable modes up to 
m=12 for Model Ml, for example, and the wavelengths of their high harmonics are found to 
be shorter than the scale height. It should be noted, however, that they have smaller growth 
rates and subdominant in driving the instability. Note also that the above analysis neither 
approves nor disproves of a particular mechanism in a mathematically rigorous sense. We 
need further investigation definitely. 



6. Implications for Astrophysical Phenomena 

In the previous sections, we have found that the standing shock wave in the accretion flow 
into the Schwarzschild black hole is generally unstable to the non-axisymmetric perturbations 
and that it oscillates with large amplitudes in the non-linear regime. Here we consider the 
astrophysical implications of the shock instability, picking up Black Hole QPOs and GRBs 
as examples. 



6.1. Black Hole QPOs 



As mentioned already, quasi-periodic oscillations have been observed for a couple of black 
hole candidates and they are attributed to some activities of the accretion disk around the 
black hol e. The shock osci lla tion model for black hole QPOs has been investigated by man y 
authors JDas et al.l l2003al B Ichakrabarti et.al.l 120041 : koki et all I2004J : bkuda et all l2007h . 
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Recently, for example, lOkuda et all (120071 ) performed two-dimensional pseudo-Newtonian 
numerical simulations of the shock oscillation in the meridian section assuming axisymmetry 
and taking into account the cooling and heating of gas and the radiation transport. They 
demonstrated that a quasi-periodically oscillating shock wave is formed around a black hole. 
They compared the numerical results with the observations for GRS 1915+105 and suggested 
that the intermedia te frequency QPO of this source might be due to the shock instability 
ffokuda et al.ll2007h . 



Our models are different from those of lOkuda et al.l (120071 ). They consider the axisym- 
metric oscillations whereas we investigate the non-axisymmetric oscillations. We take general 
relativity fully into account. Besides, we calculate the energy density spectra for the present 
models. In so doing, we employ the data in the non-linear regime, that is, 100 ms after the 
onset of computations. Note that the dynamics in the non-linear phase is almost identical in 
all the models, including the ones, in which the m = 2 or 3 mode is initially imposed instead 
of the m = 1 mode. 

Figure [TT] shows the power spectra for the m = 0, 1, 2, 3 modes in Model Ml. It is found 
that the m = mode has a quasi-periodic feature around 8 Hz, which corresponds to the 
period of large oscillations observed for the m = mode. Although there are some hints of 
other QPOs, they are much less remarkable. T his axisymmetric quasi-periodic oscillation 
is similar to those found by lOkuda et al.l (120071 ). The most important point here, however, 
is the fact that the quasi-periodicity of m = mode is induced by the non-axisymmetric 
instability through the quadratic mode coupling. 

Similar quasi-periodic oscillations are found also in other models. Their frequencies de- 
pend on the unperturbed flow and, hence, on the Bernoulli constant and specific angular mo- 
mentum. The quantitative comparison with observations is beyond the scope of this paper, 
since we have neglected radiative processes, viscosity and, among other things, we have con- 
sidered only the equatorial plane. It can be mentioned, however, that the non-axisymmetric 
shock instability is a good candidate of the source of QPOs and further investigations are 
certainly needed. 



6.2. Fluctuations in GRB jets 

Long GRBs are currently thought to be associated with massive stellar collapses and 
the subsequent formation of black holes. Although the central engine remains a mystery, it 
is widely believed that a highly relativistic jet is somehow produced near the black hole and 
its kinetic energy is later dissipated in internal shocks at larger distances, emitting gamma 
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rays (see iMeszaroa (120061 ) for a recent review). In the so-called 'patchy shell' model, the 
jet actually consists of mass shells that have slightly different velocities and collide with 
each other, generating the internal shock waves. Although the time scale of the velocity 
fluctuations is thought to be set by the dynamical time scale of the black hole, the exact 
physical processes producing the velocity variations are unknown at present. 

During the collapse of massive stars giving rise to GRBs, a large amount of matter 
will accrete on a time scale of seconds onto a proto neutron star at first and into a black 
hole later. If the progenitor is rotating rapidly prior to the collapse (jYoon et al.ll2008l ). the 
accreting matter will form a disk around th e compact object at the center. The accretion 
disk is expected to be advection-dominated (IPopham et al.lll998l ). We are thus interested in 
the stability of the accretion flows into the black hole, especially in the accretion-dominated 
regime. Here we consider the accretion flows with a standing shock wave in them, since the 
core bounce produces a shock wave, which becomes an standing shock in the core soon after 
and will continue to exist in the subsequent accretions onto a proto neutron star, the phase 
just preceding the black hole formation. Even if the bounce shock does not survive, there 
will be a lot of chances of shock formation as long as the standing shock is robust, since the 
velocity and pressure of accreting matter are fluctuating in reality. 

According to the patchy shell model, gamma rays are emitted when the kinetic energy 
of ultra-relativistic jet is dissipated in internal shock waves, which are originated from the 
inhomogeneity of the jet. Although the mechanism of jet formation remains unknown, the 
black hole is supposed to be involved. The source of the inhomogeneity is also an unsolved 
problem. If a standing shock wave exists in the accretion flow, for example, as a relic of the 
shock wave produced at the core bounce, we speculate that the intrinsic instability of the 
system against non-radial perturbations will be a natural source of fluctuations in the GRBs 
jet if it is forme d from some interact i ons b etween the accretion disk and black hole, which 
is not unlikely (IBlandford fc Znajekl 119771 ) . It is mentioned incidentally that the recent 



proge nitor models, which could produce GRB (IMacFadyen &: Woosleyl Il999l ; iHeger et al. 



20051 ). predict the injection p arameters that are appropriate for the existence of a standing 
accretion shock. For example, IHeger et al.l (120051 ) calculated the evolutions of massive stars, 
taking into account magnetic fields and obtained the specific angular momentum of several 
xlO 16 cm 2 /s and the temperature < lO 10 ^ for the matter that will later form an accretion 
disk. These numbers are just suitable for the existence of a standing shock wave in the 
accretion disk around a black hole of several M . 

Owing to the non-axisymmetric instability, the mass flux fluctuates very much indeed 
in our models. In this context, it is interesting to mention that the quasi-periodic oscilla- 
tion with a period much longer than the dynamical time scale, which we have found in the 
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previous sections, may leave its imprint somehow in the prompt gamma ray emissions or 
earl y X-ray afterglows. I t is certainly necessary, however, to study possible effects of cool- 
ing (IPopham et al.lll998l ) on the instability. The disk thickness we ignored in this paper is 
also a concern in the future work. We finally mention that the gravitational radiation by 
the non-radial shock instability may also have interesting implications. 



7. Summary and Conclusion 

We have investigated the non-axisymmetric shock instability in the accretion disk around 
Schwarzschild black holes, employing the fully general relativistic hydrodynamic simulations 
as well as the linear analysis. Both the linear and non-linear phases have been analyzed in 
detail. We have also given some possible implications for astrophysically interesting phe- 
nomena such as Black Hole QPOs and GRBs. 

The main findings in the present work £1X6 cLS follows: 

(a) The standing shock is generally unstable against non-axisymmetric perturbations, 
and a spiral arm structure is formed as a result of the growth of instability. It is typically 
one-armed, implying that the dominant mode in the non-linear phase is the m = 1 mode. 

(b) In the linear phase, the dynamical simulations are in good agreement with the linear 
analysis in such features as stability, oscillation and growth time scales. The progressive 
modes, in which the deformed shock pattern rotates in the same direction as the unperturbed 
flow, are unstable and the retrogressive modes are stable. This is consistent with the previous 
works. 

(c) In the non-linear phase, various modes are produced by non-linear couplings, which 
are mainly of quadratic nature, and the amplitudes are saturated. The axisymmetric mode 
is also induced by the non-axisymmetric instability, and the shock radius oscillates with 
large amplitudes. The oscillation periods become slightly longer than in the linear analysis 
because of larger shock radii. 

(d) Even though strong perturbations are added initially, the shock remains to exist. 
Hence the disk plus shock system is quite robust in this sense. 

(e) The comparison of various cycle time scales with the linear growth times seems to 
support the claim that the instability is induced by the acoustic- acoustic cycle, although the 
inner reflection point is not identified unambiguously. It is important to note in this respect 
that the wavelength of perturbations is longer than the scale height, which does not allow 
the WKB approximation. 
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(f) The Black Hole SASI found by lMolteni et al.l (119991 ) may be a promising candidate 
for the sources of the Black Hole QPOs and fluctuations in GRB jets. 

In the present study, we have also found that the non-axisymmetric instability is sensi- 
tive to the structure of the unperturbed steady flow. The general relativity is important in 
this respect. It should be stressed that the injection parameters that allow the existence of 
a standing shock wave are different between the GR and pseudo-Newtonian treatments. In 
fact, we have found by the direct comparison that the maximum specific angular momentum 
for the existence of multiple sonic points is different by more than 60% for the Bernoulli 
constant E < 1.003. 

Note also that the general relativity is indispensable in discussing the accretion into 
a Kerr blac k hole, since the frame-dragg ing will play an important role. This is currently 
undertaken (INagakura k Yamada 1 120081 ) . For more detailed comparison with observations, 
it is necessary to include the cooling and heating for GRBs case, and the magnetic field and 
viscosity for Black Hole QPOs. Last but not least, the discussed simulations including the 
polar dimension are inevitable. 



We are grateful to Kenta Kiuchi and Yu Yamamoto for useful discussions. This work 
was partially supported by the Grant-in-Aid for the 21st century COE program "Holistic 
Research and Education Center for Physics of Self-organizing Systems" of Waseda University 
and for Scientific Research of the Ministry of Education, Science, Sports and Culture of Japan 
(17540267, 14079202). 



A. General Relativistic Hydrodynamic Code 



Here we describe the GRHD code that are used in this paper. As mentioned already, 
it is base on the so-called central scheme, which guarantees a good accuracy even if flows 
i nclude strong shocks and /or high Lorentz factors. Magnetic fields can b e also included 
JPelZanna k BucciantiAooi Ishibata k Sekiguchill2005l : buez et al-lboosh . 



Though we do not take into account the evolution of gra vitational fie l d, th e so-called 
3+1 formalism is suitable for hydrodynamics as well. Following iDuez et al.l (120051 ). we write 
the metric in the form 



ds 2 = -a 2 dt 2 + 7ij (dx i + ftdi) (dx j + ftdt) 



(Al) 



where a, j3\ and 7^ are the lapse, shift vector, and spatial metric, respectively. The basic 
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equations for fluid dynamics in 3+1 form are expressed as: 

d tP * + dj (p*v j ) = 0, (A2) 

dtSi + d^ay/jTl) = l -a^T^g a(3)l) (A3) 

d t T + d i (a?SyT**-p.v i ) = s, (A4) 

where various variables are defined as follows: 

v* = ~ v (A5) 
u 

P* = a^/l Pou\ (A6) 

Sj = a-y/7 T° = p*hiij, (A7) 

r = afi/jT 00 - p* = p*aW - ^p - p*, (A8) 

s = a^7 { (T 00 /3 4 /? j + 2T 04 /? J ' + r j ) ^ - (T m l3 l + T 0i ) . (A9) 

In the above equations, 7 and ify are the determinant of the three metric and extrinsic 
curvature, respectively. We refer to p*, Sj and r as "conserved variables (collectively denoted 
by U)" , whereas po, p and v l are called "primitive variables (collectively expressed as P)". 

The conserved variables can be calculated directly from the primitive variables via 
Eqs. flA6j) . (IA7I) and (IA8I) . There is no analytical expression for the primitive variables 
as a function of the conserved variables, on the other hand. Since we update the conserved 
variables rather than the primitive variables, we must need to solve the latter numerically 
at each time step because they are necessary for the calculations of the characteristic wave 
speed at each cell interfa ce as shown later. If we use a T-law EOS, the inversion can be con- 
ducted easily as done by lDuez et al.l (120051 ). The same method can not be applicable to the 
general EOS, however. Hence we take a different procedure based on the Newton-Raphson 
method, which will be explained below. 

We first write down a useful relation between u l and Uj 

. 1 
u l = -jl + 7 ii w i w i } 2 . (A10) 

We define two more quantities as 



h = po^p^+^SiSjj-p^h 2 , (All) 
f 2 = r + p* - p^ahu 1 + a/7P- (A12) 



We then search iteratively for the primitive variables that satisfy /1 = f 2 — 0. We first 
guess two thermodynamical quantities po and p. Then other thermodynamical quantities 
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can be obtained from the EOS. Next, we obtain Uj from Eq. flATj) using Sj, p* and h, u f 
is determined by Uj from Eq. flAlOp . Thus the right hand sides of Eqs. flAlip and (IA12I) 
are expressed as a function of only two thermodynamical quantities. We solve them by the 
Newton-Raphson method. The initial guess is obtained from the values at the previous step. 

The net flux at the cell interface is given by the approximate solution to the Riemann 
problem. Our code adopts the HLL (Harten. Lax, and van Leer) flux, which does not require 
the complete knowledge of the solutions to the Riemann problem but the maximum wave 
speed in each direction is needed. The first step for calculating the flux is to obtain P R and 
Pl-i which are the values of primitive variables interpolated to the ri ght- and left-hand side 



of each cell interface. We have implemented the MUSCL method fjHIRSCH. C.I Il990l ) for 



this purpose. From P R and Pl, the maximum wave spe ed on each side of the cell interface, 



c± ) r and c-t L, can be calculated as in lDuez et al.l (120051 ) 



The HLL flux is then expressed with the maximum wave speeds defined by c +max = 
max(0, c +iR , c +tL ) and C- max = max(0, c_ R , c_ L ) as 

r _ ^—maxfK C+maxflj C—max^+max (Ur ^ l) I \ 

Jint — — j l Alc V 

C-mai ~r C+max 

where f R and fi are the fluxes calculated with P R and Pl, respectively. Note that if we define 
c_ max = c +max = max(0,c +ijR ,c +jL ,c_ jR ,c_ iL ), then f int becomes the local Lax-Friedrichs 
flux. 
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Table 1. Model Parameters 





A di cibcit i c 


Bernoulli 


Specific Angular 


Inner Sonic 


Shock Point 


Mach Number 


Initial 


Initicil Pcrturbeition 


IVIOUCI 


Index r 




Momentum A 


Point ri nso 


r sh 




Perturbation Mode 


Amplitude 


Ml 


4 


1.004 




!j . <J 1 V 1 3|t 


16 1 

J. \j . J. i v i * 


2.4 


I 


1 % 


M2 


I 








OQ o A/f 






1 % 


M3 


1.004 


3.50M, 


5.0M* 


34.8M, 


2.1 


1 


1 % 


M4 


1 


1.004 




4.8M* 


78 liW. 


1.5 




1 % 


M5 


I 


1.001 


3.50A4* 


5.1M* 


16.9M* 


4.1 


1 


1 % 


M6 


3 


1.005 


3.50M* 


5.0M* 


50.2M, 


1.6 


1 


1 % 


M7 


1.033 


1.13 


3.80M, 


4.4M, 


38.7M, 


2.2 


1 


1 % 


M8 


1.167 


1.02 


3.70M, 


4.6M, 


64.2M, 


1.4 


1 


1 % 


M9 


1.167 


1.02 


3.60M* 


5.0M, 


14.0M, 


2.7 


1 


1 % 


M10 


1.167 


1.03 


3.60M* 


5.0M* 


32.4M* 


1.5 


1 


1 % 


Mil 


1.433 


1.001 


3.35M* 


5.2M* 


40.6M* 


2.3 


1 


1 % 


M12 


1.433 


1.004 


3.15M* 


6.0M, 


36.5M, 


1.3 


1 


1 % 


Mlm2 


4 

1 

1 


1.004 


3.43A4* 


5.3M, 


16.1M, 


2.4 


2 


1 % 


Mlm3 


1.004 


3.43M, 


5.3M* 


16.1M» 


2.4 


O 


1 % 


MlalO 




1.004 


3.43M* 


5.3M* 


16.1M* 


2.4 


1 


10 % 


MlalOO 


3 


1.004 


3.43M* 


5.3M* 


16.1M* 


2.4 


1 


100 % 



Note. — The locations of inner sonic point and shock surface are determined by the adiabatic index, Bernoulli constant and specific angular momentum. 
The mach number is calculated in the corotating observer's frame. M, is the mass of the central black hole. 



Table 2. Properties of Instability 





dominant mode 


maximum amplitude 






Model 


in non-linear phase 


of m = mode 


Tla 




Ml 


m=l 


3.9 


Si 100ms 


ss 20ms 


M2 


m=l 


5.1 


ss 120ms 


«20ms 


M3 


m=l 


3.2 


Si200ms 


ss20ms 


M4 


m=l or 2 


1.3 




ss 50ms 


M5 


m=l 








M6 


m=l 


1.5 


ss210ms 


Si30ms 


M7 


m=l 


4.5 




ss50ms 


M8 


m=l 


1.4 


Si300ms 


Si20ms 


M9 


m=l 


3.5 


ss80ms ss20ms 




M10 


m=l 


1.3 




SilOms 


Mil 


m=l 


4.2 


«300ms 


ss60ms 


M12 










Mlm2 


m=l 


3.1 


as 100ms 


«20ms 


Mlm3 


m=l 


2.9 


Si 100ms 


«20ms 


MlalO 


m=l 


3.5 


Si 100ms 


Sd20ms 


MlalOO 


m=l 


4.0 


Si 100ms 


«20ms 



Note. — ria (ism) is the large- (small-) amplitude oscillation period. The symbol 
(-) implies no identifications. 



Table 3. Cycle Frequencies 





Corotation Point 


Oscillation Period 


Growth Time 


Wavelength of Acoustic Perturbations 










Model 


7°coro 


tosci 




Ait? 


T d~c(c) 


^c — c(c) 


r d-c(s) 


^c — c(s) 


Ml 


47.7km (10.6M*) 


3.7ms 


3.2ms 


142.6km (31.7M*) 


3.0ms 


2.0ms 


6.9ms 


5.1ms 


M2 


64.3km (14.3M*) 


6.4ms 


4.8ms 


213.7km (47. 5M*) 


5.0ms 


3.2ms 


10.6ms 


7.1ms 


M3 


62.3km (13.9A2*) 


5.9ms 


8.2ms 


171.9km (38.2Al»J 


12.yms 


7.8ms 


18.3ms 


11.3ms 


M4 


61.2km (13.6Af») 


5.6ms 


29.0ms 


120.1km (26.7M.) 


49.7ms 


31.1ms 


55.1ms 


33.9ms 


M5 


50.4km (11.2Af») 


4.1ms 


2.4ms 


147.6km (32.8M*) 


3.3ms 


1.9ms 


7.5ms 


4.9ms 


M6 


53.1km (11.8M*) 


4.5ms 


17.3ms 


114.3km (25.4Af*) 


25.0ms 


16.1ms 


29.0ms 


18.6ms 


M7 


104.8km (23. 3M*) 


14.5ms 


8.4ms 


294.3km (65.4A/*) 


20.8ms 


7.9ms 


159.2ms 


18.1 ms 


M8 


61.2km (13.6M*) 


5.5ms 


22.6ms 


121.5km (27.0M*) 


46.0ms 


26.1ms 


56.6ms 


30.0ms 


M9 


43.6km (9.7M») 


3.1ms 


2.2ms 


97.2km (21.6M*) 


3.1ms 


1.8ms 


7.5ms 


4.9ms 


M10 


61.2km (13. 6M*) 


5.6ms 


11.0ms 


157.9km (35.1M.) 


12.8ms 


8.5ms 


19.3ms 


12.3ms 


Mil 


67.5km (15.0M.) 


7.2ms 


10.7ms 


211.0km (46. 9M.) 


14.6ms 


9.3ms 


19.6ms 


12.8ms 


M12 


44.1km (9.8Af.) 


3.5ms 


252.2ms 


113.8km (25.3JW.) 


16.1ms 


13.0ms 


18.9ms 


15.6ms 


Mlm2 


51.7km (11.5M*) 


2.2ms 


2.7ms 


82.8km (18.4Af«) 


2.5ms 


1.6ms 


6.9ms 


5.1ms 


Mlm3 


54.0km (12.0M.) 


1.5ms 


2.6ms 


58.9km (13.1M*) 


2.2ms 


1.5ms 


6.9ms 


5.1ms 



-J 

Note. — tosaii tgrow/2'K and \ w represent the oscillation period, growth time and wavelength of acoustic perturbations, respectively, which are obtained I 
by linear analysis. r [J _ c ( c ), t c _ c ( c ), t < i_ c i s ), t c _ c ( s ) are obtained from Eqs. H43H to H46I I, which show the acoustic-acoustic cycle or advective-acoustic 
cycle between the shock surface and the corotation or inner sonic point. 
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radius[Mj I [Mj 



Fig. 1. — Left: The location of sonic points as a function of the Bernoulli constant (E) and 
specific angular momentum (A) in the Schwarzschild geometry. The solid, dashed, dotted 
lines correspond to E — 1.004, 1.02 and 1.1, respectively. Right: The injection parameters 
for the existence of a standing shock wave. The shaded region allows the standing shock 
wave. The adiabatic index is 4/3 for both panels. 
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t = (ms) t = 21.5(ms) t = 42.0(ms) 



t = 70.0(ms) t = HO.O(ms) t = 262.0(ms) 




t = 318.0(ms) t = 425.0(ms) 



Fig. 2. — The time evolution of velocity for Model Ml. The color contour shows the mag- 
nitude of radial velocity. The arrows represent the velocities at their positions. The central 
region in blue corresponds to the black hole. 
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I I I I 

t = (ms) t = 21.5(ms) t = 42.0(ms) 

rm 

t = 70.0(ms) t = HO.O(ms) t = 262.0(ms) 




t = 318.0(ms) t = 425.0(ms) 



Fig. 3. — The time evolution of entropy for Model Ml. 




t [msj t[ms] 

Fig. 4. — The time evolutions of the m = mode (left) and the m — 1,2,3 modes (right). 
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4,0e+02 




0,0e+00 



O.Oe+00 5.0e+03 1.0e+04 1 .5e+04 2.O&+04 
w r [l/sj 



Fig. 5. — The real and imaginary parts of eigen frequencies for some of the m = 1, 2, 3 modes 
for Model Ml. 




Fig. 6. — The comparison of the time evolutions of the amplitudes of m = 1 mode for 
Model Ml obtained by the linear analysis and dynamical simulation. The red line shows the 
evolution expected by the linear analysis, while the green crosses are the simulation results. 
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4,0e+02 




O.Oe+00 



O.Oe+00 4.0e+03 8.Oe+03 Ue+04 1.6e+04 
w r [1/s] 



4,0e+02 




O.Oe+00 

1.0e+03 



1.2e+03 1.4e+03 1.6e+03 
w r /m [1/s] 



1.8e+03 



Fig. 8. — Left: the most unstable eigenfrequency in each mode. Right: the same as in left 
figure but the horizontal axis is the pattern frequency u r /m 
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Fig. 9. — The time evolutions of the amplitudes of various modes. The upper left (right) 
panel shows the results for Model Mlm2 (Mlm3). The lower panel displays the evolution of 
m = mode for Models Ml, Mlm2 and Mlm3. 



3 
5,5 

2 
15 

1 

5 





ad-actcor) 
ac-ac[corJ 



J* 
f 



4 — 

3.5 - 

3 - 

2.S - 



0.1 




ad-ac(so) 
ac-ac(so} 



Fig. 10. — The ratio of the growth rate to the frequencies of advective-acoustic (+) and 
acoustic-acoustic cycles (x) for all the models. In the left (right) panel, the corotation 
(inner sonic) point is assumed to be the inner reflection point. 
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Fig. 11. — Power spectra of energy density for the m = mode (left) and m — 1, 2, 3 modes 
(right). 



